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PREFACE 

The  work  reported  herein  was  conducted  by  the  Arnold  Engineering  Development 
Center  (AEDC),  Air  Force  Systems  Command  (AFSC).  The  Air  Force  project  manager  was 
Elton  R.  Thompson,  DOT.  The  results  of  the  research  were  obtained  by  ARO,  Inc.,  AEDC 
Group  (a  Sverdrup  Corporation  Company),  operating  contractor  for  the  AEDC,  AFSC, 
Arnold  Air  Force  Station,  Tennessee,  under  ARO  Project  No.  P32A-01.  The  manuscript 
was  submitted  for  publication  on  September  16,  1980. 

The  authors  acknowledge  the  support  and  contributions  of  W.  E.  Dietz,  T.  L.  Donegan, 
R.  R.  Jones,  K.  R.  Stansbury,  T.  W.  Swafford,  and  D.  L.  Whitfield  in  demonstrating  the 
accuracy  and  flexibility  of  the  computer  program. 
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1.0  INTRODUCTION 

Recent  success  has  been  demonstrated  in  the  difficult  area  of  computing  transonic,  three- 
dimensional  inviscid  flow  fields  over  complex  configurations.  Notable  is  the  work  of 
Jameson  and  Caughey  (Refs.  1  through  3)  for  solution  of  the  full  potential  equations  with 
wing-body  configurations,  and  Boppe  (Ref.  4)  for  solution  of  the  small-perturbation- 
potential  equations  with  wing-body-pylon-nacelle  configurations.  Another  major 
contribution  is  evident  in  the  work  of  Cline  (Ref.  5),  who  solved  the  axisymmetric  Euler 
equations  for  convergent-divergent  nozzle  configurations.  Although  as  yet  unpublished, 
Phares  of  ARO,  Inc.,  has  extended  the  Cline  code  to  three-dimensional  nozzle 
configurations. 

These  computer  programs  are  used  routinely  at  AEDC  to  resolve  problems  and  to 
enhance  the  ground  test  capabilities  of  the  wind  tunnel  and  engine  test  facilities.  Because 
many  situations  arise  for  which  the  available  codes  are  not  applicable,  it  was  deemed 
necessary  to  write  a  new  computer  program  with  as  few  limitations  as  possible.  Specifically, 
there  were  to  be  no  geometry  restrictions,  and  the  code  should  be  capable  of  treating  both 
internal  and  external  configurations  without  the  assumptions  of  irrotational  or  isoenergetic 
flow. 

This  report  describes  the  development  of  a  computer  program  (designated  ARO-1)  to 
solve  the  three-dimensional,  unsteady  Euler  equations  in  Cartesian  coordinates  using  a  finite 
volume  (volume  flux)  approach.  The  basic  numerical  algorithm  is  the  explicit  predictor- 
corrector  scheme  of  MacCormack  (Ref.  6).  Some  results  are  compared  with  exact  solutions 
to  illustrate  the  basic  accuracy  of  the  code,  whereas  other  solutions  are  presented  to 
demonstrate  the  great  flexibility  available  with  respect  to  geometry  and  extreme  flow 
situations. 


2.0  NUMERICAL  METHOD 


2.1  EULER  EQUATIONS 

The  three-dimensional,  unsteady  Euler  equations  are  conservatively  written  in  Cartesian 
coordinates  as 


<9G 

dt 


+  V  •  F 


0 


(1) 
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and  the  pressure  is  given  by 

p  =  (y-i: 


—  p  ^u2  +  v2  +  w2  ^ 


(3) 


All  variables  are  dimensionless,  with  the  reference  conditions  usually  taken  to  be  free-stream 
density  and  sound  speed,  the  latter  given  by 


•-(Vp) 

Integrating  Eq.  (1)  over  a  small  stationary  volume,  V,  yields 


—  f  GdV  +  f 
9t  J  V  Jv 


V  •  F  dV  =  0 


Using  the  mean  value  and  divergence  theorems  results  in 


dG 

dt 


(4) 


(5) 


(6) 


where  n  is  the  outward  unit  normal  to  the  surface,  S,  enclosing  V,  and  G  is  associated  with 
some  point  interior  to  the  volume  (taken  to  be  the  volume  center  for  practical  purposes). 


2.2  COMPUTATIONAL  MESH 

The  grid  imposed  on  the  computational  flow  region  is  topologically  equivalent  to  a 
regular  grid  on  a  cube.  The  resulting  control  volumes  are  hexahedrons  with  quadrilateral 
faces  which  are  not  necessarily  planar.  Degenerate  volumes  may  also  be  used;  for  example, 
two  edges  may  coincide  to  form  a  distorted  triangular  prism,  thereby  maintaining  flexibility 
with  respect  to  arbitrary  geometry.  A  representative  control  volume  within  the  mesh  is 
shown  in  Fig.  1.  Three  pairs  of  opposite  faces  are  evident  and  denoted  by  subscripts  k  =  1, 
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2,  or  3.  One  of  the  eight  vertices  is  designated  as  the  positive  corner,  and  the  area  vectors  of 
the  three  faces  meeting  at  the  vertex  are  designated  bySi+,  S^,  and  S/.  The  other  three 
faces  have  area  vectors  designated  by  Sf,  Sf,  and  S3-,  respectively.  When  two  volume 
elements  share  a  common  face,  that  face  is  represented  by  in  one  volume  and  in  the 
other.  This  coupling  forms  three  “pseudo  directions”  through  the  computational  mesh  and 
permits  selection  of  three  indices  for  cataloging  the  flow  variable  independent  of  the  mesh 
orientation  with  respect  to  the  coordinate  axes. 


The  area  vector  is  computed  as  one-half  the  cross  product  of  the  diagonal  vectors.  The 
resulting  direction  is  an  “average”  normal  over  the  surface  with  the  magnitude  being  the 
projected  area  in  that  direction.  This  vector  is  resolved  into  Cartesian  components.  Center 
points  are  defined  by  simple  averaging;  the  center  point  coordinates  of  each  face  are  the 
average  of  the  four  vertex  coordinates,  and  the  coordinate  of  the  volume  center  is  the 
average  of  the  eight  vertex  coordinates.  A  center  vector  is  defined  as  the  vector  directed  from 
the  volume  center  to  the  center  of  a  face.  The  volume  of  each  hexahedron  is  calculated  as 
one-third  the  sum  of  the  dot  products  of  the  center  and  area  vectors. 

Routines  for  mesh  generation  are  user-supplied.  Several  examples  of  usable 
computational  meshes  are  given  in  Section  3.0.  The  only  restriction  with  respect  to  mesh 
coordinates  is  that  they  be  consistently  ordered  so  that  the  resulting  volume  calculation  is 
positive.  Negative  volumes  indicate  either  errors  in  constructing  the  mesh  or  simply  a 
reversal  of  one  of  the  coordinate  indices. 

2.3  SOLUTION  ALGORITHM 

The  basic  numerical  algorithm  is  the  explicit  predictor-corrector  method  of 
MacCormack  (Ref.  6).  The  flux  integral  in  Eq.  (6)  is  approximated  by  summing  the  scalar 
products  of  the  area  vectors  and  the  appropriately  evaluated  flux  vectors.  The  solution  can 
be  advanced  from  time  n  to  n  -I- 1  using 


—  n  +  1 

G 


n 

G 


Cl)  A  t 
V 


F‘st  + 


(7) 
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where  F  is  the  flux  vector  evaluated  at  the  volume  center,  F+  is  the  flux  vector  evaluated  at 
the  center  point  of  the  volume  sharing  the  face,  and  F~  similarly  corresponds  with  ST- 
The  overbar  represents  evaluation  at  the  predictor  time.  Pressure  is  updated  using  Eq.  (3) 
and  boundary  conditions  applied  at  the  end  of  both  predictor  and  corrector  steps.  The 
orientation  of  the  predictor-corrector  is  changed  each  time  step  by  cycling  the  positive 
corner  among  the  eight  vertices. 

The  parameter,  co,  was  introduced  in  Eq.  (7)  to  help  stabilize  the  algorithm.  For  to  =  1, 
the  scheme  reduces  to  the  volume  flux  form  of  the  MacCormack  algorithm.  Increasing  co 
increases  the  truncation  error,  which  acts  as  an  artificial  viscosity,  thereby  smoothing  the 
results.  Increasing  <o  also  rigorously  reduces  the  order  of  accuracy  from  second  order  to  first 
order.  However,  for  an  irregular  computational  mesh  the  algorithm  is  less  than  second-order 
accurate  even  with  co  =  1 . 


The  time  step,  At,  is  limited  in  magnitude  by  the  Courant-Friedrich-Lewy  stability 
criterion.  The  maximum  allowable  time  step  is  computed  as  the  minimum  of 


At 


_ V 

I  q  •  s  |  +  c  |  s 


(8) 


over  the  entire  mesh  for  each  face  of  each  volume  where  q  is  the  velocity  vector.  Preliminary 
experience  indicates  that  the  time  step  must  be  further  reduced  by  for  stability.  Since  the 
sound  speed,  c,  requires  a  square-root  operation,  At  is  calculated  only  every  128th  time  step 
to  reduce  computer  time. 


2.4  BOUNDARY  AND  INITIAL  CONDITIONS 


Boundary  conditions  are  applied  by  using  phantom  points  exterior  from  the 
computational  mesh.  When  the  properties  such  as  upstream  free-stream  conditions  are 
known,  these  values  are  specified  at  the  phantom  points. 


Mirror  conditions  are  used  at  planes  of  symmetry.  Pressure,  density,  and  energy  at  the 
phantom  points  are  equated  to  the  neighboring  interior  values.  The  velocity  vector  is 
mirrored  using 


q 


be 


q 

in 


(9) 


where  q;n  is  the  velocity  vector  neighboring  the  phantom  point,  qt>c>  and  S  is  the  boundary 
area  vector. 
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At  wall  or  body  surface  boundaries,  mirror  conditions  are  also  generally  used.  This 
results  in  inaccuracies  when  the  boundaries  are  not  planar;  however,  no  efficient  method  of 
implementing  a  normal  momentum  equation  or  characteristic-type  boundary  condition  has 
been  developed  for  arbitrary  surface  geometry.  The  magnitude  of  the  errors  introduced  by 
the  mirror  assumption  is  discussed  in  Section  3.0. 

Outflow  or  outer  boundaries  are  treated  with  either  the  first  or  second  differences  of  all 
variables  set  to  zero.  In  general,  zero  gradient  conditions  obtained  by  equating  the  phantom 
points  to  the  neighboring  interior  values  have  proved  to  be  more  stable  than  zero  second 
difference  conditions. 

Initial  conditions  for  exterior-type  flows  are  usually  taken  to  be  free-stream  values  with 
body  surface  boundary  conditions  impulsively  applied  at  time  zero.  For  interior  flows  the 
initial  conditions  are  usually  obtained  from  one-dimensional  isentropic  flow  relations  with 
the  velocity  vector  aligned  with  the  mesh.  Subsonic  duct  flow  solutions  also  use  an  upstream 
boundary  condition  which  is  similar  to  that  of  Cline  (Ref.  5),  derived  from  characterisitc 
relations  with  pressure  specified  at  the  outflow  boundary. 

2.5  SMOOTHING 

During  development  of  the  code  it  became  very  evident  that  the  basic  algorithm  is 
unconditionally  unstable  when  applied  to  an  irregular  mesh.  Some  form  of  artificial 
viscosity  is  required  to  damp  the  high  frequency  truncation  errors  and  to  dissipate  expansion 
shocks  which  are  admissible  (Ref.  7).  Three  types  of  smoothing  have  been  used. 

The  “omega  factor”  (co)  introduced  in  Eq.  (7)  is  a  viable  method  of  favorably  increasing 
the  truncation  error  to  add  artificial  viscosity.  Its  primary  advantage  is  the  addition  of 
stability  with  minimal  impact  on  efficiency.  However,  the  omega  factor  has  been  rarely  used 
during  the  code  evaluation  for  two  reasons:  (1)  aesthetically,  it  seems  incorrect  to  reduce  the 
formal  accuracy  of  a  numerical  scheme,  and  more  to  the  point,  (2)  the  production  version  of 
the  computer  program  contained  a  coding  error  relative  to  the  omega  factor. 

Some  solutions  have  been  obtained  using  the  third-order  smoothing  of  Ref.  8.  This 
procedure  specifically  evaluates  the  third-order  truncation  errors  and  includes  them  in 
evaluation  of  the  surface  flux  integral  over  each  volume  element.  Although  stability  is 
assured  and  the  solutions  are  accurate,  the  computation  time  is  significantly  increased 
relative  to  the  basic  algorithm. 
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The  majority  of  solutions  have  been  obtained  using  explicit  weighted  averaging.  At  each 
eighth  time  step,  the  entire  flow  field  is  smoothed  by  simple  averaging  of  adjacent  volume 
elements,  with  the  central  element  being  weighted  by  the  factor 

W  =  2° 

\  min 


where  Vmin  represents  the  minimum  of  volume  V  over  the  computational  domain.  After 
averaging,  the  smoothed  interior  values  are  used  to  update  boundary  conditions. 

2.6  CODE  STRUCTURE  AND  TIMING  STUDIES 

The  basic  code  consists  of  six  subroutines  which  accomplish  the  following  tasks: 

1.  DTCAL  computes  the  maximum  allowable  time  step. 

2.  BC  applies  boundary  conditions. 

3.  FLUXX  evaluates  surface  fluxes  in  the  k  =  1  direction. 

4.  FLUXY  evaluates  surfaces  fluxes  in  the  k  =  2  direction. 

5.  FLUXZ  evaluates  surface  fluxes  in  the  k  =  3  direction. 

6.  UPDATE  advances  the  solution  in  time. 

For  efficiency  of  computation,  Eq.  (7)  is  coded  in  the  following  form 


—  n  +  1 

G 

n 

G  - 

(l>  At/ 

v  V 

£  F  •  S  + 
k  k 

+  ?  F'  ■  ST 

k  K 

n  +  l 

—  n  +  1 

At 

I7r  - + 

G 

G 

Lf 

*  St  +  L  F  • 

2V 

!_\ k 

k  k 

-  (2cu 

-nil 

-  F  •  S  + 

Lr  K 

+  £  F '  •  S  -k 
w  k 

(11) 


so  that  instead  of  storing  two  time  values  of  each  variable,  one  time  value  and  a  flux 
accumulator  are  stored  for  each  of  the  five  dependent  variables. 


The  flux  calculations  of  the  form  F»S  occur  in  pairs  in  Eq.  (11)  such  that  each  of  the 
FLUX(k)  routines  could  be  optimized  in  the  k  direction.  For  example,  the  F»S^  term  for 
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one  volume  is  the  same  as  the  F~*ST  term  for  the  adjacent  volume.  When  these  terms  are 
summed  in  the  k  direction,  unnecessary  duplication  of  computations  are  avoided.  This 
process  is  identical  to  MacCormack’s  split  operator  (Ref.  9),  except  that  all  fluxes  are 
accumulated  prior  to  updating. 

The  split  coding  necessitated  that  a  pressure  array  be  saved  to  avoid  undesirable 
repetition  of  computations.  For  the  same  reason  the  volume  and  area  vector  arrays  are  also 
saved.  Thus,  21  variables  are  stored  for  each  mesh  point.  This  number  could  be  reduced  by 
six  if  the  area  vectors  from  the  mesh  coordinates  were  recalculated  as  required,  but  the 
computation  time  would  increase  by  about  25  percent. 

Most  of  the  results  presented  were  obtained  from  the  Cray-1  computer  at  United 
Computer  Services,  Kansas  City.  The  Cray  operating  system  contains  a  flow  trace  routine 
which  monitors  the  time  required  in  each  subroutine.  The  resulting  (mesh-dependent) 
timings  are  given  in  the  following  table: 

* 

Routine  Percent  of  Time 


DTCAL 

6 

BC 

17 

FLUXX 

16 

FLUXY 

13 

FLUXZ 

19 

UPDATE 

14 

Other 

15 

A  count  of  the  number  of  operations  per  time  step  (predictor  plus  corrector)  within 
FLUX(k)  and  UPDATE  shows  112  floating  point  additions  and  108  floating  point 
multiplications  per  mesh  point.  Central  processor  time  averages  about  7  x  10~6  sec  per  mesh 
point  per  time  step,  depending  on  the  mesh;  this  means  the  Cray-1  is  operating  at  50  million 
floating  point  operations  per  second  when  it  is  in  these  four  routines. 

3.0  RESULTS  AND  DISCUSSION 

To  ensure  the  accuracy  and  proper  coding  of  the  computer  program,  the  flow  field  about 
several  simple  configurations  was  computed  for  comparison  with  results  which  have  exact 
solutions.  The  geometries  included  a  two-dimensional  double-wedge  and  a  cone  at  angle  of 
attack  at  supersonic  conditions.  Two-dimensional  planar  (or  axisymmetric)  configurations 
are  represented  using  one  planar  (or  wedge-shaped)  array  of  volume  elements  with  symmetry 
boundary  conditions  on  two  sides. 
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Solutions  were  obtained  for  Mach  number  2.0  flow  over  a  4.5-deg,  double-symmetric 
wedge  using  two  computational  grids.  A  basically  orthogonal  mesh  yields  the  results  given  in 
Fig.  2.  Surface  pressures  are  in  good  agreement  with  theory,  but  results  within  the  flow  field 
are  not.  Dispersion  of  the  shock  waves  is  evident,  particularly  so  in  the  contour  plot  of 
constant  pressure  lines  of  Fig.  2b.  Results  obtained  with  the  mesh  aligned  with  the  leading- 
and  trailing-edge  shocks  are  given  in  Fig.  3.  The  surface  pressures  are  in  excellent  agreement 
with  theory,  and  little  dispersion  or  dissipation  is  evident.  However,  mesh  alignment  permits 
the  appearance  of  an  expansion  shock  emanating  from  the  apex  of  the  wedge.  This 
phenomenon  is  attributable  to  the  availability  of  an  expansion  shock  as  an  exact  solution  to 
the  Euler  equations  across  one  mesh  boundary  at  the  expansion  location  (Ref.  7).  Some 
form  of  dissipation  is  necessary  to  prevent  the  occurrence  of  expansion  shocks.  Results 
obtained  with  third-order  and  explicit  smoothing  are  presented  in  Fig.  4.  It  should  be  noted 
that  the  “omega  factor”  is  nondissipative  and  thus  has  little  effect  on  the  expansion  shock 
for  this  mesh. 

Computations  were  also  made  for  Mach  number  2.0  flow  over  a  40-deg  cone  at  5-deg 
angle  of  attack.  The  computational  mesh  is  illustrated  in  Fig.  5a,  and  an  isobar  contour  plot 
is  given  in  Fig.  5b.  Some  smearing  of  the  bow  shock  is  apparent.  Stability  problems  were 
encountered  with  the  cone  mesh  configuration  so  that  smoothing  was  mandatory.  Surface 
pressure  results  are  given  in  Fig.  5c  for  both  the  mesh  of  Fig.  5a  and  for  a  coarser  mesh 
containing  only  nine  circumferential  elements.  Both  solutions  are  in  good  agreement  with 
the  results  of  Jones  (Ref.  10).  However,  detailed  examination  of  the  solutions  shows  the 
fine-mesh  solution  with  less  than  twice  the  number  of  points  to  be  about  five  times  more 
accurate  than  the  coarse  mesh  solution. 

It  is  relatively  easy  to  incorporate  streamline  curvature  effects  in  the  surface  boundary 
condition  for  simple  configurations.  Computations  have  been  made  for  a  1.5-  caliber  ogive- 
cylinder  at  20-deg  angle  of  attack  with  both  mirror  and  the  normal  momentum  equation 
surface  boundary  conditions.  The  results  are  given  in  Fig.  6,  along  with  experimental  data 
from  Ref.  11.  Cross  flow  separation  occurs  at  x/D  »  4.5  (off  the  figure  scale)  which,  of 
course,  cannot  be  predicted  using  the  Euler  equations.  Before  separation,  the  computations 
and  experimental  data  are  in  reasonably  good  agreement. 

Although  only  speculation  at  present,  the  error  introduced  through  use  of  the  mirror 
boundary  condition  appears  to  create  effective  vorticity  and  a  corresponding  stagnation 
pressure  loss  that  bears  some  resemblance  to  true  viscous  effects.  For  example,  calculations 
indicate  a  5-percent  loss  in  stagnation  pressure  close  to  the  ogive-cylinder  body  at  x/D  =  2 
which  increases  to  10-percent  loss  at  x/D  =  10.  The  loss  is  circumferentially  uniform  close 
to  the  nose  but  then  convected  to  the  leeward  side,  as  expected  of  viscous  effects.  Users  of 
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the  code  should  routinely  compute  and  examine  both  stagnation  pressure  and  enthalpy 
distribution  throughout  the  flow  field  to  evaluate  solution  accuracy. 

An  example  of  the  program’s  flexibility  is  indicated  by  the  mesh  in  Fig.  7a  for  computing 
the  flow  over  three  separated  aircraft  stores  in  a  triple-ejection-rack  configuration.  The 
mesh  consists  of  two  distinct  computational  grids  interconnected  by  the  boundary  condition 
routine,  which  required  the  addition  of  only  35  FORTRAN  statements  to  the  basic  computer 
program.  A  comparison  of  computed  results  with  the  experimental  surface  pressure  data  of 
Heltsley  and  Cline  (Ref.  12)  is  given  in  Fig.  7b.  Note  that  the  flow  in  the  gap  between  the 
stores  is  locally  supersonic.  Computations  at  a  free-stream  Mach  number  of  0.9  (not  shown) 
indicate  two  strong  shocks  in  the  gap. 

As  a  precursor  to  the  development  of  a  three-dimensional  wing/cascade  version  of 
ARO-1,  solutions  have  been  obtained  for  a  two-dimensional  lifting  airfoil.  The  calculations 
are  compared  with  the  experimental  data  of  Ref.  11  in  Fig.  8.  The  trailing-edge  boundary 
condition  was  simply  continuity  of  dependent  variables  across. the  “wake.”  The  lack  of 
agreement  between  the  computations  and  experiment,  particularly  at  supercritical 
conditions,  is  unexplained.  However,  these  differences  are  quite  similar  to  the  effects  of 
wind  tunnel  wall  interference  illustrated  in  Ref.  11. 

Computation  of  flow  over  an  F-16  aircraft  forebody  (including  strakes  but,  as  yet, 
without  inlet  flow)  provides  a  three-dimensional  example.  The  computational  mesh  in  the 
vertical  plane  is  given  in  Fig.  9a,  and  a  comparison  of  computed  results  with  experimental 
pressure  data  along  the  top  surface  is  presented  in  Fig.  9b.  Considering  the  relative 
coarseness  of  the  mesh,  these  results  clearly  show  the  general  applicability  of  ARO-1  to  this 
class  of  configurations. 

Use  of  the  Euler  equations  rather  than  the  potential  flow  assumption  permits  calculation 
of  rotational  flows  such  as  secondary  flow  induced  in  curved  ducts.  Such  a  calculation  can 
be  applied  in  the  contraction  section  of  a  wind  tunnel.  The  outermost  part  of  the 
computational  mesh  for  a  contraction  with  transition  from  a  circular  to  a  rectangular  to  a 
square  cross  section  is  illustrated  in  Fig.  10a,  with  the  interior  mesh  at  the  upstream  cross 
section  given  in  Fig.  10b.  With  uniform  flow  at  the  inlet,  the  exit  flow  is  contaminated  by 
small  secondary  flows,  as  shown  by  the  cross  flow  velocity  vectors  in  Fig.  10c.  The 
maximum  magnitude  of  the  corresponding  flow  angularity  is  0.1  deg. 

Another  example  of  a  curved  duct  computation  is  given  in  Fig.  11.  An- aircraft  engine- 
inlet  duct  is  modeled  with  transition  from  an  elliptic  to  a  circular  cross  section  with  a  mild 
centerline  offset.  Note  the  replicated  mesh  points  at  the  duct  centerline  in  Fig.  lib  and  recall 
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that  the  finite  volume  formulation  does  not  permit  flux  through  a  surface  of  zero  area. 
Nonetheless,  computed  results  such  as  those  given  in  Fig.  lie  appear  reasonable  and  are  in 
good  qualitative  agreement  with  experimental  data. 

A  final  example  of  code  flexibility  is  provided  by  solutions  of  an  axisymmetric  flow  over 
a  hollow-nose  infrared  sensing  probe  as  shown  in  Fig.  12a.  The  mesh  (Fig.  12b)  was 
generated  by  a  variant  of  Thomas’  approach  (Ref.  13)  which  solves  two  coupled  Poisson 
equations  to  yield  interior  mesh  coordinates  from  specified  boundary  coordinates. 
Computed  isobars  for  steady  flow  at  M*  =  1.35  are  given  in  Fig.  12c.  The  explicit 
smoothing  added  to  maintain  stability  tends  to  smear  the  bow  shock  over  several  mesh 
points.  Cancellation  of  the  shock  at  the  mesh  boundary  was  obtained  using  zero  normal 
gradient  boundary  conditions. 

Time-accurate  unsteady  flow  can  also  be  calculated,  as  evidenced  by  the  sequence  of 
isobar  plots  in  Fig.  13  resulting  from  the  introduction  of  a  planar  blast  wave  upstream  of  the 
infrared  sensing  probe.  Shock-shock  interactions  and  multiple  shock  reflections  pose  no 
computational  difficulties,  although  the  accuracy  of  the  computations  is  unknown,  except 
to  note  that  the  blast  wave  speed  was  correctly  represented. 

The  probe  was  to  be  cooled  by  injecting  air  in  the  cavity.  A  solution  was  obtained  with, 
in  effect,  counter-flowing  cold  jets  by  local  modification  of  the  surface  boundary  condition 
to  permit  blowing.  Isobars  of  the  resulting  solution  are  given  in  Fig.  14a,  and  the  path  of  the 
“coolant”  is  indicated  by  lines  of  constant  temperature  in  Fig.  14b.  The  face  of  the  probe 
cannot  be  cooled  by  this  technique. 

4.0  CONCLUDING  REMARKS 

A  generalized  computer  program  for  solution  of  the  three-dimensional,  unsteady  Euler 
equations  has  been  developed.  Selection  of  the  finite-volume  approach  in  Cartesian 
coordinates  yielded  a  flexible  code  capable  of  treating  arbitrary  model  geometry,  including 
both  internal  and  external  flow  configurations.  Efficient  coding  and  use  of  the  Cray-1 
computer  enables  reasonably  short  solution  times  on  fairly  complex  and  dense  meshes. 

Development  of  the  code  is  providing  even  greater  flexibility  and  extending  the  number 
of  model  flow  configurations  that  can  be  computed  on  a  routine  basis.  For  example, 
viscous-inviscid  interactions  can  be  computed  approximately  through  incorporation  of  an 
inverse,  two-dimensional  strip,  boundary-layer  method  (Ref.  14)  which  requires 
modification  of  only  the  boundary  condition  routine.  That  code  is  designated  ARO-2. 
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a.  Computational  mesh 


b.  Isobars 

Figure  2.  Computation  of  flow  over  a  double-wedge  with  orthogonal  mesh. 
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d.  Pressure  at  y  =  0,2 
Figure  2.  Concluded, 
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Figure  3.  Computation  of  flow  over  a  double-wedge  with  inclined  mesh. 
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Figure  5.  Computation  of  flow  over  a  cone  at  angle  of  attack. 
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c.  Surface  pressure 
Figure  5.  Concluded. 
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a.  Computational  mesh 

Figure  7.  Computation  of  flow  over  a  three-store  cluster 
at  zero  angle  of  attack. 
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a.  Computational  mesh 

Figure  8.  Computation  of  flow  over  a  lifting  airfoil. 
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x/c 

d.  Surface  pressure,  =  0.73 
Figure  8.  Concluded. 
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c.  Exit  plane  velocity  vectors 

Figure  10.  Computation  of  flow  within  a  wind  tunnel  contraction. 
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b.  Computational  mesh 

Figure  12.  Computation  of  flow  about  a  hollow-nose  probe. 
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c.  Isobars 

Figure  12.  Concluded. 
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e.  Time  =  5  At 


f.  Time  =  6  At 


Figure  13.  Computation  of  unsteady  flow. 
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NOMENCLATURE 

c  Sound  speed 

D  Diameter 

e  Specific  energy 

F  Tensor  of  independent  variables  defined  in  Eq.  (2) 

G  Vector  of  dependent  variables,  Eq.  (2) 

k  Direction  index 

n  Unit-normal  vector 

P  Pressure 

q  Velocity  vector 

5  Surface 

Sk  Surface  area  vector 

t  Time 

u  Velocity  in  x-direction 

V  Volume 

v  Velocity  in  y-direction 

W  Weighting  factor,  Eq.  (10) 

w  Velocity  in  z-direction 

x,y,z  Cartesian  coordinates 

a  Angle  of  attack 

7  Ratio  of  specific  heats 

6  Density 

Roll  orientation 

(a  Stabilization  parameter  in  Eq.  (7) 
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SUBSCRIPTS 

be  Boundary  condition 

in  Interior 

k  Direction  index 

max  Maximum 

min  Minimum 

oo  Free-stream  conditions 

SUPERSCRIPTS 
n  Time  index 

+  Positive  direction  sense 

Negative  direction  sense 
overbar  Predictor  time  level 
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